
In this exercise you will learn and understand
Convolution is a form of simulating multispectral data from hyperspectral data via the Relative Spectral Responses (RSR). RSRs define the shape of the multispectral sensors and resulting data.


For example to generate the NIR band (above) we integrate all the hyperspectral data in the range defined by the RSR of the NIR band as would the actual multispectral sensor do.
The spectral convolution process uses the following equation:
Where:
The spatial convolution process uses the following equation:
Where:
The spatial convolution is an aggregation process of the finer resolution pixles to the coarser ones. For example to create a 250m resolution pixel from 1 m pixels we simply need to average 250x250 pixels window. In practice there are other ways to doing this, that we will cover later.
# load library modules
import os
import numpy as np
import matplotlib.pyplot as plt
import viplab_lib3 as vip
# convolution library
import viplab_convolution as vipConv
import h5py
import io
from PIL import Image
# This code forces plotting inline within the code (below the cell) and not outsoide in separate windows
%matplotlib inline
import warnings
warnings.filterwarnings('ignore')
conda install -c conda-forge pyhdf
Collecting package metadata (current_repodata.json): done
Solving environment: done
==> WARNING: A newer version of conda exists. <==
current version: 4.11.0
latest version: 4.12.0
Please update conda by running
$ conda update -n base -c defaults conda
## Package Plan ##
environment location: /Users/condongroup/opt/miniconda3/envs/remsens
added / updated specs:
- pyhdf
The following packages will be downloaded:
package | build
---------------------------|-----------------
certifi-2021.10.8 | py38h50d1736_2 145 KB conda-forge
openssl-1.1.1n | h6c3fc93_0 1.9 MB conda-forge
pyhdf-0.10.3 | py38he785675_1 145 KB conda-forge
python_abi-3.8 | 2_cp38 4 KB conda-forge
------------------------------------------------------------
Total: 2.2 MB
The following NEW packages will be INSTALLED:
hdf4 conda-forge/osx-64::hdf4-4.2.15-hefd3b78_3
pyhdf conda-forge/osx-64::pyhdf-0.10.3-py38he785675_1
python_abi conda-forge/osx-64::python_abi-3.8-2_cp38
The following packages will be SUPERSEDED by a higher-priority channel:
ca-certificates pkgs/main::ca-certificates-2022.3.29-~ --> conda-forge::ca-certificates-2021.10.8-h033912b_0
certifi pkgs/main::certifi-2021.10.8-py38hecd~ --> conda-forge::certifi-2021.10.8-py38h50d1736_2
openssl pkgs/main::openssl-1.1.1n-hca72f7f_0 --> conda-forge::openssl-1.1.1n-h6c3fc93_0
Downloading and Extracting Packages
openssl-1.1.1n | 1.9 MB | ##################################### | 100%
python_abi-3.8 | 4 KB | ##################################### | 100%
certifi-2021.10.8 | 145 KB | ##################################### | 100%
pyhdf-0.10.3 | 145 KB | ##################################### | 100%
Preparing transaction: done
Verifying transaction: done
Executing transaction: done
Note: you may need to restart the kernel to use updated packages.
hdf5_file = 'NEON_D14_SRER_DP3_502000_3523000_reflectance.h5'
hdf5= h5py.File(hdf5_file, "r")
#list_dataset lists the names of datasets in an hdf5 file
def list_dataset(name,node):
if isinstance(node, h5py.Dataset):
print(name)
hdf5.visititems(list_dataset)
SRER/Reflectance/Metadata/Ancillary_Imagery/Aerosol_Optical_Depth SRER/Reflectance/Metadata/Ancillary_Imagery/Aspect SRER/Reflectance/Metadata/Ancillary_Imagery/Cast_Shadow SRER/Reflectance/Metadata/Ancillary_Imagery/Dark_Dense_Vegetation_Classification SRER/Reflectance/Metadata/Ancillary_Imagery/Data_Selection_Index SRER/Reflectance/Metadata/Ancillary_Imagery/Haze_Cloud_Water_Map SRER/Reflectance/Metadata/Ancillary_Imagery/Illumination_Factor SRER/Reflectance/Metadata/Ancillary_Imagery/Path_Length SRER/Reflectance/Metadata/Ancillary_Imagery/Sky_View_Factor SRER/Reflectance/Metadata/Ancillary_Imagery/Slope SRER/Reflectance/Metadata/Ancillary_Imagery/Smooth_Surface_Elevation SRER/Reflectance/Metadata/Ancillary_Imagery/Visibility_Index_Map SRER/Reflectance/Metadata/Ancillary_Imagery/Water_Vapor_Column SRER/Reflectance/Metadata/Ancillary_Imagery/Weather_Quality_Indicator SRER/Reflectance/Metadata/Coordinate_System/Coordinate_System_String SRER/Reflectance/Metadata/Coordinate_System/EPSG Code SRER/Reflectance/Metadata/Coordinate_System/Map_Info SRER/Reflectance/Metadata/Coordinate_System/Proj4 SRER/Reflectance/Metadata/Logs/175032/ATCOR_Input_file SRER/Reflectance/Metadata/Logs/175032/ATCOR_Processing_Log SRER/Reflectance/Metadata/Logs/175032/Shadow_Processing_Log SRER/Reflectance/Metadata/Logs/175032/Skyview_Processing_Log SRER/Reflectance/Metadata/Logs/175032/Solar_Azimuth_Angle SRER/Reflectance/Metadata/Logs/175032/Solar_Zenith_Angle SRER/Reflectance/Metadata/Logs/175509/ATCOR_Input_file SRER/Reflectance/Metadata/Logs/175509/ATCOR_Processing_Log SRER/Reflectance/Metadata/Logs/175509/Shadow_Processing_Log SRER/Reflectance/Metadata/Logs/175509/Skyview_Processing_Log SRER/Reflectance/Metadata/Logs/175509/Solar_Azimuth_Angle SRER/Reflectance/Metadata/Logs/175509/Solar_Zenith_Angle SRER/Reflectance/Metadata/Logs/175923/ATCOR_Input_file SRER/Reflectance/Metadata/Logs/175923/ATCOR_Processing_Log SRER/Reflectance/Metadata/Logs/175923/Shadow_Processing_Log SRER/Reflectance/Metadata/Logs/175923/Skyview_Processing_Log SRER/Reflectance/Metadata/Logs/175923/Solar_Azimuth_Angle SRER/Reflectance/Metadata/Logs/175923/Solar_Zenith_Angle SRER/Reflectance/Metadata/Logs/180334/ATCOR_Input_file SRER/Reflectance/Metadata/Logs/180334/ATCOR_Processing_Log SRER/Reflectance/Metadata/Logs/180334/Shadow_Processing_Log SRER/Reflectance/Metadata/Logs/180334/Skyview_Processing_Log SRER/Reflectance/Metadata/Logs/180334/Solar_Azimuth_Angle SRER/Reflectance/Metadata/Logs/180334/Solar_Zenith_Angle SRER/Reflectance/Metadata/Spectral_Data/FWHM SRER/Reflectance/Metadata/Spectral_Data/Wavelength SRER/Reflectance/Metadata/to-sensor_azimuth_angle SRER/Reflectance/Metadata/to-sensor_zenith_angle SRER/Reflectance/Reflectance_Data
#ls_dataset displays the name, shape, and type of datasets in hdf5 file
def ls_dataset(name,node):
if isinstance(node, h5py.Dataset):
print(node)
# Use the visititems methods to learn abotu the data
hdf5.visititems(ls_dataset)
# Look at the the last line printed, that is teh data we will work with
<HDF5 dataset "Aerosol_Optical_Depth": shape (1000, 1000), type "<i2"> <HDF5 dataset "Aspect": shape (1000, 1000), type "<f4"> <HDF5 dataset "Cast_Shadow": shape (1000, 1000), type "|u1"> <HDF5 dataset "Dark_Dense_Vegetation_Classification": shape (1000, 1000), type "|u1"> <HDF5 dataset "Data_Selection_Index": shape (1000, 1000), type "<i4"> <HDF5 dataset "Haze_Cloud_Water_Map": shape (1000, 1000), type "|u1"> <HDF5 dataset "Illumination_Factor": shape (1000, 1000), type "|u1"> <HDF5 dataset "Path_Length": shape (1000, 1000), type "<f4"> <HDF5 dataset "Sky_View_Factor": shape (1000, 1000), type "|u1"> <HDF5 dataset "Slope": shape (1000, 1000), type "<f4"> <HDF5 dataset "Smooth_Surface_Elevation": shape (1000, 1000), type "<f4"> <HDF5 dataset "Visibility_Index_Map": shape (1000, 1000), type "|u1"> <HDF5 dataset "Water_Vapor_Column": shape (1000, 1000), type "<f4"> <HDF5 dataset "Weather_Quality_Indicator": shape (1000, 1000, 3), type "|u1"> <HDF5 dataset "Coordinate_System_String": shape (), type "|O"> <HDF5 dataset "EPSG Code": shape (), type "|O"> <HDF5 dataset "Map_Info": shape (), type "|O"> <HDF5 dataset "Proj4": shape (), type "|O"> <HDF5 dataset "ATCOR_Input_file": shape (), type "|O"> <HDF5 dataset "ATCOR_Processing_Log": shape (), type "|O"> <HDF5 dataset "Shadow_Processing_Log": shape (), type "|O"> <HDF5 dataset "Skyview_Processing_Log": shape (), type "|O"> <HDF5 dataset "Solar_Azimuth_Angle": shape (), type "<f4"> <HDF5 dataset "Solar_Zenith_Angle": shape (), type "<f4"> <HDF5 dataset "ATCOR_Input_file": shape (), type "|O"> <HDF5 dataset "ATCOR_Processing_Log": shape (), type "|O"> <HDF5 dataset "Shadow_Processing_Log": shape (), type "|O"> <HDF5 dataset "Skyview_Processing_Log": shape (), type "|O"> <HDF5 dataset "Solar_Azimuth_Angle": shape (), type "<f4"> <HDF5 dataset "Solar_Zenith_Angle": shape (), type "<f4"> <HDF5 dataset "ATCOR_Input_file": shape (), type "|O"> <HDF5 dataset "ATCOR_Processing_Log": shape (), type "|O"> <HDF5 dataset "Shadow_Processing_Log": shape (), type "|O"> <HDF5 dataset "Skyview_Processing_Log": shape (), type "|O"> <HDF5 dataset "Solar_Azimuth_Angle": shape (), type "<f4"> <HDF5 dataset "Solar_Zenith_Angle": shape (), type "<f4"> <HDF5 dataset "ATCOR_Input_file": shape (), type "|O"> <HDF5 dataset "ATCOR_Processing_Log": shape (), type "|O"> <HDF5 dataset "Shadow_Processing_Log": shape (), type "|O"> <HDF5 dataset "Skyview_Processing_Log": shape (), type "|O"> <HDF5 dataset "Solar_Azimuth_Angle": shape (), type "<f4"> <HDF5 dataset "Solar_Zenith_Angle": shape (), type "<f4"> <HDF5 dataset "FWHM": shape (426,), type "<f4"> <HDF5 dataset "Wavelength": shape (426,), type "<f4"> <HDF5 dataset "to-sensor_azimuth_angle": shape (1000, 1000), type "<f4"> <HDF5 dataset "to-sensor_zenith_angle": shape (1000, 1000), type "<f4"> <HDF5 dataset "Reflectance_Data": shape (1000, 1000, 426), type "<i2">
hdf5
<HDF5 file "NEON_D14_SRER_DP3_502000_3523000_reflectance.h5" (mode r)>
surf_refl = hdf5['SRER']['Reflectance']
print(surf_refl)
<HDF5 group "/SRER/Reflectance" (2 members)>
reflArray = surf_refl['Reflectance_Data']
print(reflArray)
## Notice the data dimension and shape
<HDF5 dataset "Reflectance_Data": shape (1000, 1000, 426), type "<i2">
refl_shape = reflArray.shape
print('Reflectance Data Dimensions:',refl_shape)
Reflectance Data Dimensions: (1000, 1000, 426)
mapInfo = surf_refl['Metadata']['Coordinate_System']['Map_Info']
print('Data Cube Map Info: ',mapInfo[()])
#First convert mapInfo to a string
mapInfo_string = str(mapInfo[()]) #convert to string
#split the strings using the separator comma ","
mapInfo_split = mapInfo_string.split(",")
print(mapInfo_split)
#Extract the upper left-hand corner coordinates from mapInfo
xMin = float(mapInfo_split[3])
yMax = float(mapInfo_split[4])
#Calculate the xMax and yMin values from the dimensions using the array size and pixel size
xMax = xMin + (refl_shape[1]*res[0]) #xMax = left edge + (# of columns * x pixel resolution)
yMin = yMax - (refl_shape[0]*res[1]) #yMin = top edge - (# of rows * y pixel resolution)
#Define extent as a tuple:
cube_ext = (xMin, xMax, yMin, yMax)
print('cube_ext: ',cube_ext)
print('cube_ext type: ',type(cube_ext))
Data Cube Map Info: b'UTM, 1.000, 1.000, 502000.00, 3524000.0, 1.0000000, 1.0000000, 12, North, WGS-84, units=Meters, 0' ["b'UTM", ' 1.000', ' 1.000', ' 502000.00', ' 3524000.0', ' 1.0000000', ' 1.0000000', ' 12', ' North', ' WGS-84', ' units=Meters', " 0'"]
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Input In [24], in <module> 10 yMax = float(mapInfo_split[4]) 12 #Calculate the xMax and yMin values from the dimensions using the array size and pixel size ---> 13 xMax = xMin + (refl_shape[1]*res[0]) #xMax = left edge + (# of columns * x pixel resolution) 14 yMin = yMax - (refl_shape[0]*res[1]) #yMin = top edge - (# of rows * y pixel resolution) 16 #Define extent as a tuple: NameError: name 'res' is not defined
</center>
#define the wavelengths variable
wavelengths = surf_refl['Metadata']['Spectral_Data']['Wavelength']
#View wavelength information and values
print('wavelengths:',wavelengths)
wavelengths: <HDF5 dataset "Wavelength": shape (426,), type "<f4">
# Display min & max wavelengths
print('min wavelength:', np.amin(wavelengths),'nm')
print('max wavelength:', np.amax(wavelengths),'nm')
## Notice that we are not having to define these manually and that we are using what we call metdata to get this info
min wavelength: 381.5437 nm max wavelength: 2509.932 nm
#show the band widths between the first 2 bands and last 2 bands
print('band width between first 2 bands =',(wavelengths[1]-wavelengths[0]),'nm')
print('band width between last 2 bands =',(wavelengths[-1]-wavelengths[-2]),'nm')
band width between first 2 bands = 5.0079956 nm band width between last 2 bands = 5.0080566 nm
mapInfo = surf_refl['Metadata']['Coordinate_System']['Map_Info']
print('Data Cube Map Info:',mapInfo[()])
Data Cube Map Info: b'UTM, 1.000, 1.000, 502000.00, 3524000.0, 1.0000000, 1.0000000, 12, North, WGS-84, units=Meters, 0'
#First convert mapInfo to a string
mapInfo_string = str(mapInfo[()]) #convert to string
#split the strings using the separator ","
mapInfo_split = mapInfo_string.split(",")
print(mapInfo_split)
["b'UTM", ' 1.000', ' 1.000', ' 502000.00', ' 3524000.0', ' 1.0000000', ' 1.0000000', ' 12', ' North', ' WGS-84', ' units=Meters', " 0'"]
#Extract the resolution & convert to floating decimal number
res = float(mapInfo_split[5]),float(mapInfo_split[6])
print('Resolution:',res)
Resolution: (1.0, 1.0)
#Extract the upper left-hand corner coordinates from mapInfo
xMin = float(mapInfo_split[3])
yMax = float(mapInfo_split[4])
#Calculate the xMax and yMin values from the dimensions
xMax = xMin + (refl_shape[1]*res[0]) #xMax = left edge + (# of columns * x pixel resolution)
yMin = yMax - (refl_shape[0]*res[1]) #yMin = top edge - (# of rows * y pixel resolution)
#Define extent as a tuple:
cube_ext = (xMin, xMax, yMin, yMax)
print('cube_ext:',cube_ext)
print('cube_ext type:',type(cube_ext))
cube_ext: (502000.0, 503000.0, 3523000.0, 3524000.0) cube_ext type: <class 'tuple'>
Band_56 = reflArray[:,:,55].astype(float)
print('b56 type:',type(Band_56))
print('b56 shape:',Band_56.shape)
print('Band 56 Reflectance:\n',Band_56)
b56 type: <class 'numpy.ndarray'> b56 shape: (1000, 1000) Band 56 Reflectance: [[ 720. 787. 1125. ... 1580. 983. 840.] [ 929. 1289. 1289. ... 1661. 989. 989.] [ 765. 1038. 1475. ... 1928. 1080. 670.] ... [1383. 1425. 1757. ... 3972. 4721. 4826.] [1344. 1475. 1711. ... 4904. 4711. 4631.] [1565. 1709. 924. ... 4288. 4378. 3744.]]
#View and apply scale factor and data ignore value
scaleFactor = reflArray.attrs['Scale_Factor']
noDataValue = reflArray.attrs['Data_Ignore_Value']
print('Scale Factor:',scaleFactor)
print('Data Ignore Value:',noDataValue)
Band_56[Band_56==int(noDataValue)]=np.nan
Band_56 = Band_56/scaleFactor
print('Cleaned Band 56 Reflectance:\n',Band_56)
Scale Factor: 10000.0 Data Ignore Value: -9999.0 Cleaned Band 56 Reflectance: [[0.072 0.0787 0.1125 ... 0.158 0.0983 0.084 ] [0.0929 0.1289 0.1289 ... 0.1661 0.0989 0.0989] [0.0765 0.1038 0.1475 ... 0.1928 0.108 0.067 ] ... [0.1383 0.1425 0.1757 ... 0.3972 0.4721 0.4826] [0.1344 0.1475 0.1711 ... 0.4904 0.4711 0.4631] [0.1565 0.1709 0.0924 ... 0.4288 0.4378 0.3744]]
data_plot = plt.imshow(Band_56,extent=cube_ext,cmap='Greys')
plt.hist(Band_56[~np.isnan(Band_56)],50); #50 signifies the # of bins
### We can see that most of the reflectance values are within the range 0 - 0.4
#### In order to show more contrast in the image, we can adjust the colorlimit (clim) to 0-0.4 (or 0 to 0.2 only)
data_plot = plt.imshow(Band_56,extent=cube_ext,cmap='gist_earth',clim=(0,0.4))
plt.title('Data Cube Band 56 Reflectance');
# Assign some common bands by name and number. Feel free to experiment
# and use other regions of the spectrum that correspond to the color of interest.
# Review the references slides in the Lab. presentation about where these band numbers come from
bandRED=48
bandGREEN=34
bandBLUE=17
bandNIR=97
# Combine the Red, Green and Blue data into an RGB image for display
print("Creating RGB Image...")
RGBImage=vip.Image_getRGB(reflArray[:,:,bandRED],reflArray[:,:,bandGREEN],reflArray[:,:,bandBLUE],5000)
# Display RGB True color Image
plt.figure(figsize=(10,10))
plt.title('RGB True color')
plt.axis('off')
plt.imshow(RGBImage)
Creating RGB Image...
<matplotlib.image.AxesImage at 0x7f9865778670>
# create Spectral Convolution object
NEON=vipConv.SpecConvolution()
#Load NEON wavelength values that match each NEON band
#WavesList=NEON.Load_Wavesfile('./Data/NEON_wavelength_values.txt')
# Will use the list loaded above
WavesList=wavelengths
#Load simulated sensor RSR function
NEON.Load_RSR('NEON_Sentinel2A_RSR.csv')
#display RSR basic information
NEON.displayRSR()
#Retrieve RSR for a specific band
S2_RSR_RED=NEON.getBand_RSR('RED')
S2_RSR_NIR=NEON.getBand_RSR('NIR')
#create a plot of the RSR
plt.figure()
plt.title('RSR Sensor - Put NAME Here')
plt.xlabel('Wavelength (nm)')
plt.ylabel('RSR')
#RSR_RED[:,0] is the X Axis
#RSR_RED[:,1] is the Y Axis
plt.plot(S2_RSR_RED[:,0],S2_RSR_RED[:,1],color='red')
plt.plot(S2_RSR_NIR[:,0],S2_RSR_NIR[:,1],color='black')
plt.grid(True)
plt.tight_layout()
plt.show()
Reading RSR: NEON_Sentinel2A_RSR.csv Convolution: Sensor: SENTINEL2A Resolution: 10.0 meters Bands: ['BLUE', 'GREEN', 'RED', 'NIR']
#apply convolution to dataset
SENTINELCube=NEON.Spectral_Spatial_Convolution(WavesList,reflArray)
Processing convolution... band: BLUE min= 440.0 max= 535.0 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 . band: GREEN min= 537.0 max= 582.0 32 33 34 35 36 37 38 39 40 41 . band: RED min= 646.0 max= 684.0 53 54 55 56 57 58 59 60 61 . band: NIR min= 773.0 max= 908.0 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 . Spatial Resampling at 10.0 [ 100 , 100 ]
RGBImage=vip.Image_getRGB(SENTINELCube[:,:,2],SENTINELCube[:,:,1],SENTINELCube[:,:,0],8000)
# Display RGB True color Image
plt.figure(figsize=(10,10))
plt.title('SENTINEL simulated RGB at actual Sentinel spatial resolution')
plt.imshow(RGBImage)
<matplotlib.image.AxesImage at 0x7f984bce9d00>
# Rerun Convolution at a different spatial resolution, 5 meters
SENTINELCube=NEON.Spectral_Spatial_Convolution(WavesList,reflArray,5)
RGBImage=vip.Image_getRGB(SENTINELCube[:,:,2],SENTINELCube[:,:,1],SENTINELCube[:,:,0],8000)
Processing convolution... band: BLUE min= 440.0 max= 535.0 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 . band: GREEN min= 537.0 max= 582.0 32 33 34 35 36 37 38 39 40 41 . band: RED min= 646.0 max= 684.0 53 54 55 56 57 58 59 60 61 . band: NIR min= 773.0 max= 908.0 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 . Spatial Resampling at 5 [ 200 , 200 ]
# Display RGB True color Image
plt.figure(figsize=(10,10))
plt.title('SENTINEL simulated RGB at hypothetical 5 meters')
plt.imshow(RGBImage)
<matplotlib.image.AxesImage at 0x7f984bcf4b80>
# ReRun only spectral Convolution (not spatial)
SENTINELCube=NEON.Spectral_Convolution(WavesList,reflArray)
RGBImage=vip.Image_getRGB(SENTINELCube[:,:,2],SENTINELCube[:,:,1],SENTINELCube[:,:,0],8000)
Processing convolution... band: BLUE min= 440.0 max= 535.0 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 . band: GREEN min= 537.0 max= 582.0 32 33 34 35 36 37 38 39 40 41 . band: RED min= 646.0 max= 684.0 53 54 55 56 57 58 59 60 61 . band: NIR min= 773.0 max= 908.0 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 .
# Display RGB True color Image
plt.figure(figsize=(10,10))
plt.title('SENTINEL simulated but no spatial convolution')
plt.imshow(RGBImage,extent=cube_ext)
<matplotlib.image.AxesImage at 0x7f984d995850>
#Load simulated sensor RSR function
NEON.Load_RSR('NEON_Landsat8_RSR.csv')
#display RSR basic information
NEON.displayRSR()
#Retrieve RSR for a specific band
L8_RSR_RED=NEON.getBand_RSR('RED')
L8_RSR_NIR=NEON.getBand_RSR('NIR')
#Load simulated sensor RSR function
NEON.Load_RSR('NEON_MODIS_RSR_New.csv')
#display RSR basic information
NEON.displayRSR()
#Retrieve RSR for a specific band
M_RSR_RED=NEON.getBand_RSR('RED')
M_RSR_NIR=NEON.getBand_RSR('NIR')
#Load simulated sensor RSR function
NEON.Load_RSR('NEON_VIIRS_RSR.csv')
#display RSR basic information
NEON.displayRSR()
#Retrieve RSR for a specific band
V_RSR_RED=NEON.getBand_RSR('RED')
V_RSR_NIR=NEON.getBand_RSR('NIR')
Reading RSR: NEON_Landsat8_RSR.csv Convolution: Sensor: LANDSAT8/OLI Resolution: 30.0 meters Bands: ['BLUE', 'GREEN', 'RED', 'NIR'] Reading RSR: NEON_MODIS_RSR_New.csv Convolution: Sensor: MODIS/AQUA Resolution: 250.0 meters Bands: ['BLUE', 'GREEN', 'RED', 'NIR'] Reading RSR: NEON_VIIRS_RSR.csv Convolution: Sensor: VIIRS Resolution: 500.0 meters Bands: ['BLUE', 'GREEN', 'RED', 'NIR']
#create a plot of the RSR
plt.figure(figsize=(15,5))
plt.title('RSR Sensor - Landsat 8')
plt.xlabel('Wavelength (nm)')
plt.ylabel('RSR')
#RSR_RED[:,0] is the X Axis
#RSR_RED[:,1] is the Y Axis
plt.plot(M_RSR_RED[:,0],M_RSR_RED[:,1],'-+',color='red',label='RED MODIS')
plt.plot(M_RSR_NIR[:,0],M_RSR_NIR[:,1],'-+',color='black',label='NIR MODIS')
plt.plot(V_RSR_RED[:,0],V_RSR_RED[:,1],'-o',color='red',label='RED VIIRS')
plt.plot(V_RSR_NIR[:,0],V_RSR_NIR[:,1],'-o',color='black',label='NIR VIIRS')
plt.plot(L8_RSR_RED[:,0],L8_RSR_RED[:,1],'-.',color='red',label='RED Landsat 8')
plt.plot(L8_RSR_NIR[:,0],L8_RSR_NIR[:,1],'-.',color='black',label='NIR Landsat 8')
plt.plot(S2_RSR_RED[:,0],S2_RSR_RED[:,1],color='red', label='RED Sentinel 2')
plt.plot(S2_RSR_NIR[:,0],S2_RSR_NIR[:,1],color='black', label='NIR Sentinel 2')
plt.grid(True)
plt.legend(loc=[0,-0.5], ncol=4)
plt.tight_layout()
plt.show()
#apply convolution to dataset
LANDSATCube=NEON.Spectral_Spatial_Convolution(WavesList,reflArray)
Processing convolution... band: BLUE min= 436.0 max= 527.0 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 . band: GREEN min= 513.0 max= 600.0 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 . band: RED min= 626.0 max= 682.0 49 50 51 52 53 54 55 56 57 58 59 60 . band: NIR min= 830.0 max= 896.0 90 91 92 93 94 95 96 97 98 99 100 101 102 103 . Spatial Resampling at 30.0 [ 33 , 33 ]
RGBImage=vip.Image_getRGB(LANDSATCube[:,:,2],LANDSATCube[:,:,1],LANDSATCube[:,:,0],8000)
# Display RGB True color Image
plt.figure(figsize=(10,10))
plt.title('LANDSAT 8 simulated RGB at actual spatial resolution')
plt.imshow(RGBImage)
<matplotlib.image.AxesImage at 0x7f984cea9160>
# Rerun Convolution at a different spatial resolution, 5 meters
LANDSATCube=NEON.Spectral_Spatial_Convolution(WavesList,reflArray,5)
RGBImage=vip.Image_getRGB(LANDSATCube[:,:,2],LANDSATCube[:,:,1],LANDSATCube[:,:,0],8000)
Processing convolution... band: BLUE min= 436.0 max= 527.0 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 . band: GREEN min= 513.0 max= 600.0 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 . band: RED min= 626.0 max= 682.0 49 50 51 52 53 54 55 56 57 58 59 60 . band: NIR min= 830.0 max= 896.0 90 91 92 93 94 95 96 97 98 99 100 101 102 103 . Spatial Resampling at 5 [ 200 , 200 ]
# Display RGB True color Image
plt.figure(figsize=(10,10))
plt.title('LANDSAT 8 simulated RGB at hypothetical 5 meters')
plt.imshow(RGBImage)
<matplotlib.image.AxesImage at 0x7f9865781d30>
# ReRun only spectral Convolution (not spatial)
LANDSATCube=NEON.Spectral_Convolution(WavesList,reflArray)
RGBImage=vip.Image_getRGB(LANDSATCube[:,:,2],LANDSATCube[:,:,1],LANDSATCube[:,:,0],8000)
Processing convolution... band: BLUE min= 436.0 max= 527.0 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 . band: GREEN min= 513.0 max= 600.0 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 . band: RED min= 626.0 max= 682.0 49 50 51 52 53 54 55 56 57 58 59 60 . band: NIR min= 830.0 max= 896.0 90 91 92 93 94 95 96 97 98 99 100 101 102 103 .
# Display RGB True color Image
plt.figure(figsize=(10,10))
plt.title('LANDSAT 8 simulated but no spatial convolution')
plt.imshow(RGBImage,extent=cube_ext)
<matplotlib.image.AxesImage at 0x7f984ced5c10>
#display a message the program ended
print("program ended.")
program ended.
In this exercise you will learn how to design and use a mask to assist with data and images analysis.
A Mask is usually a thematic (discrete values) image made up of ZEROS [0] and ONES [1], or other numbers/classes.
| Original Imager | Mask 1 | Mask 2 |
|---|---|---|
![]() |
![]() |
![]() |
The mask(s) define(s) the locations and extent of each features.
To understand the behavior (ex: spectral) of remote sensing data over a particular object/feature we can run a spectral signature (or analysis) over a single or few pixels from the feature, or better yet (Sample) average the spectral signatures of many pixels corresponding to that feature, or even all of them. But how?.
Masking provides an answer, and may contain more than a single feature and result from more than one test/method. The analysis will then be applied to only the feature(s) of interest.
# load library modules
import os
import numpy as np
import matplotlib.pyplot as plt
import viplab_lib3 as vip
#Load NEON Dataset
filename="NEON_Composed_Img.bsq"
nrows=500
ncols=500
nbands = 426
datatype=np.int16
# Read the wavelength values from the textfile
file = open('NEON_wavelength_values.txt', 'r')
Xvalues= file.readlines()
nvalues=len(Xvalues)
#convert text to number (float)
for i in range(0,nvalues):
Xvalues[i]=float(Xvalues[i])
# assign some bands with an index to easy remember them
bandRED=53
bandGREEN=37
bandBLUE=18
bandNIR=98
bandSWIR1=246
bandSWIR2=340
# Read all bands into a single DataCube
DataCube=vip.BSQ_band_read(filename,-426,nrows,ncols)
# Combine the Red, Green and Blue data into an RGB model for display
print("Creating RGB Image...")
# Extract some bands from the cube
DataRED=DataCube[:,:,bandRED]
DataGREEN=DataCube[:,:,bandGREEN]
DataBLUE=DataCube[:,:,bandBLUE]
DataNIR=DataCube[:,:,bandNIR]
# Display RGB Image
RGBImage=vip.Image_getRGB(DataRED,DataGREEN,DataBLUE,8000)
# Display RGB True color Image
plt.figure(figsize=(8,8))
plt.axis('off')
plt.title("NEON RGB - This a composite image made up of pieces")
plt.imshow(RGBImage)
nbands= 426 Creating RGB Image...
<matplotlib.image.AxesImage at 0x7f9857612640>
# Compute Stat values for a RED band
RED_avg=DataRED.mean()
RED_stdev=DataRED.std()
print("Full Image RED: avg=",RED_avg,", stdev=",RED_stdev)
Full Image RED: avg= 1118.491816 , stdev= 1341.629444353776
#Use masks to compute stats about objects only
# Automatic mask creation:
# Create a mask to select white roofs (using the RED band only for example)
# An example of a more complex mask with NDVI
# NDVI = (NIR-red)/(NIR+red)
#Notice how easy it is in Python to run complex data processing, in this case creating NDVI is one line of code
NDVI_Numerator=DataNIR-DataRED
NDVI_Denominator=DataNIR+DataRED
# In order to avoid division by zero we need to be careful, the code below does some checking befoe outputing NDVI
with np.errstate(divide='ignore', invalid='ignore'):
NDVI =np.true_divide(NDVI_Numerator,NDVI_Denominator)
NDVI[NDVI == np.inf] = 0
NDVI = np.nan_to_num(NDVI)
#Find roofs by using a red band threshold mask
MaskRoof=DataRED > 5000
fig, figarr = plt.subplots(1, 2, figsize=(12,6))
fig.tight_layout()
figarr[0].set_title("Roof Mask ")
figarr[0].axis('off')
figarr[0].imshow(MaskRoof)
# Display RGB Image
RGBImage=vip.Image_getRGB(DataRED,DataGREEN,DataBLUE,8000)
# Display RGB True color Image
figarr[1].axis('off')
figarr[1].set_title("NEON RGB - This a composite image made up of pieces")
figarr[1].imshow(RGBImage)
<matplotlib.image.AxesImage at 0x7f985750d280>
#Find water bodies by using a NIR band threshold mask
MaskWater=DataNIR < 100
fig, figarr = plt.subplots(1, 2, figsize=(12,6))
fig.tight_layout()
figarr[0].set_title("Water Mask")
figarr[0].axis('off')
figarr[0].imshow(MaskWater)
# Display RGB Image
RGBImage=vip.Image_getRGB(DataRED,DataGREEN,DataBLUE,8000)
# Display RGB True color Image
figarr[1].axis('off')
figarr[1].set_title("NEON RGB - This a composite image made up of pieces")
figarr[1].imshow(RGBImage)
<matplotlib.image.AxesImage at 0x7f984e524eb0>
# Mask Pecan trees (using NDVI band)
MaskPecan= NDVI > 0.90
fig, figarr = plt.subplots(1, 2, figsize=(12,6))
fig.tight_layout()
figarr[0].set_title("Pecan Trees Mask using NDVI < 0.9 threshold")
figarr[0].axis('off')
figarr[0].imshow(MaskPecan)
# Display RGB Image
RGBImage=vip.Image_getRGB(DataRED,DataGREEN,DataBLUE,8000)
# Display RGB True color Image
figarr[1].axis('off')
figarr[1].set_title("NEON RGB - This a composite image made up of pieces")
figarr[1].imshow(RGBImage)
<matplotlib.image.AxesImage at 0x7f9856bd58e0>
# Mask Pecan trees (using NDVI band)
MaskPecan= NDVI > 0.60
fig, figarr = plt.subplots(1, 2, figsize=(12,6))
fig.tight_layout()
figarr[0].set_title("Plants Mask using NDVI > 0.6 threshold")
figarr[0].axis('off')
figarr[0].imshow(MaskPecan)
# Display RGB Image
RGBImage=vip.Image_getRGB(DataRED,DataGREEN,DataBLUE,8000)
# Display RGB True color Image
figarr[1].axis('off')
figarr[1].set_title("NEON RGB - This a composite image made up of pieces")
figarr[1].imshow(RGBImage)
<matplotlib.image.AxesImage at 0x7f9856aa4700>
# Get the average value of the RED for the houses using the mask
# This code will travel the full image and only retains and works
# with any pixel that meets the Mask conditions
Roof_RED_avg=np.mean(DataRED[ MaskRoof ==True ])
Roof_RED_stdev=np.std(DataRED[ MaskRoof ==True ])
# Get the average value of the RED for the water pixels
Water_RED_avg=np.mean(DataRED[ MaskWater ==True ])
Water_RED_stdev=np.std(DataRED[ MaskWater ==True ])
#Using a predifined mask - from D2L
# Mask values
# 1 Water
# 2 White roofs
# 3 Walnut orchard trees
# 4 Forest trees
# Read a mask from file
MaskData1= vip.BSQ_band_read('NEON_Mask.bsq',0,500,500)
MaskData2= vip.BSQ_band_read('NEON_Mask.bsq',1,500,500)
#Display the masks
fig, figarr = plt.subplots(1, 3, figsize=(12,6))
fig.tight_layout()
figarr[0].set_title("Mask from file")
figarr[0].axis('off')
figarr[0].imshow(MaskData1)
figarr[1].set_title("Mask from file")
figarr[1].axis('off')
figarr[1].imshow(MaskData2)
# Display RGB Image
RGBImage=vip.Image_getRGB(DataRED,DataGREEN,DataBLUE,8000)
# Display RGB True color Image
figarr[2].axis('off')
figarr[2].set_title("NEON RGB - This a composite image made up of pieces")
figarr[2].imshow(RGBImage)
<matplotlib.image.AxesImage at 0x7f984adf0cd0>
# Get the average value of the RED for the water pixels
MaskPecan=MaskData2 == 3
Pecan_RED_avg=np.mean(DataRED[ MaskPecan ==True ])
Pecan_RED_stdev=np.std(DataRED[ MaskPecan ==True ])
MaskForest=MaskData2 == 4
Forest_RED_avg=np.mean(DataRED[ MaskForest ==True ])
Forest_RED_stdev=np.std(DataRED[ MaskForest ==True ])
print("RED stats using masks")
print("Houses roof: avg=",Roof_RED_avg,", stdev=",Roof_RED_stdev)
print("Water: avg=",Water_RED_avg,", stdev=",Water_RED_stdev)
print("Walnut Trees: avg=",Pecan_RED_avg,", stdev=",Pecan_RED_stdev)
print("Forest Trees: avg=",Forest_RED_avg,", stdev=",Pecan_RED_stdev)
RED stats using masks Houses roof: avg= 7200.777116919706 , stdev= 1024.063550604201 Water: avg= 747.7773073666384 , stdev= 102.12372098425168 Walnut Trees: avg= 322.8010876132931 , stdev= 178.50904048131653 Forest Trees: avg= 506.33531500288996 , stdev= 178.50904048131653
#Using the predifined mask -
# Mask values
# 1 Water
# 2 White roofs
# 3 Walnut trees
# 4 Forest
# 5 Grass
# 6 Road
# 0 everything else
Objects_Label = ['Pecan Trees','Water']
Objects_Color = ['darkgreen','blue']
Avg_Spectra=np.zeros((2,426))
for band in range(0,nbands) :
Data2D=DataCube[:,:,band]
Avg_Spectra[0,band]=np.average(Data2D[ MaskData2 ==3 ])
for band in range(0,nbands) :
Data2D=DataCube[:,:,band]
Avg_Spectra[1,band]=np.average(Data2D[ MaskData2 ==1 ])
plt.figure(figsize=(12,6))
plt.title('Features Spectral Signatures',fontsize=20)
plt.xlabel('Wavelength (nm)',fontsize=12)
plt.ylabel('Reflectance ',fontsize=12)
for i,type in enumerate(Objects_Label):
Spectra=Avg_Spectra[i]
plt.plot(Xvalues,Spectra,color=Objects_Color[i],label=Objects_Label[i])
plt.grid()
plt.legend(loc='upper left')
<matplotlib.legend.Legend at 0x7f9856936430>
# Data SWIR1 and SWIR2 have not been created yet
DataSWIR1=DataCube[:,:,bandSWIR1]
DataSWIR2=DataCube[:,:,bandSWIR2]
bands = {'Blue':DataBLUE, 'Green':DataGREEN, 'Red':DataRED,
'NIR':DataNIR, 'SWIR1':DataSWIR1, 'SWIR2':DataSWIR2}
for x,y in bands.items():
Data = y
name = x
Roof_avg=np.mean(Data[ MaskRoof ==True ])
Roof_stdev=np.std(Data[ MaskRoof ==True ])
Water_avg=np.mean(Data[ MaskWater ==True ])
Water_stdev=np.std(Data[ MaskWater ==True ])
MaskPecan=MaskData2 == 3
Pecan_avg=np.mean(Data[ MaskPecan ==True ])
Pecan_stdev=np.std(Data[ MaskPecan ==True ])
MaskForest=MaskData2 == 4
Forest_avg=np.mean(Data[ MaskForest ==True ])
Forest_stdev=np.std(Data[ MaskForest ==True ])
print(name, "band stats using masks")
print("Houses roof: avg=",Roof_avg,", stdev=",Roof_stdev)
print("Water: avg=",Water_avg,", stdev=",Water_stdev)
print("Pecan Trees: avg=",Pecan_avg,", stdev=",Pecan_stdev)
print("Forest Trees: avg=",Forest_avg,", stdev=",Pecan_stdev)
print()
Blue band stats using masks Houses roof: avg= 6334.894819220535 , stdev= 1194.442016196077 Water: avg= 374.25486875529214 , stdev= 114.5380342543871 Pecan Trees: avg= 229.69679758308158 , stdev= 76.5304814696349 Forest Trees: avg= 472.00379819998346 , stdev= 76.5304814696349 Green band stats using masks Houses roof: avg= 6853.501956487713 , stdev= 1074.4543962789367 Water: avg= 735.9729043183743 , stdev= 132.59962507414457 Pecan Trees: avg= 547.9279758308157 , stdev= 178.34344309420035 Forest Trees: avg= 697.447774750227 , stdev= 178.34344309420035 Red band stats using masks Houses roof: avg= 7200.777116919706 , stdev= 1024.063550604201 Water: avg= 747.7773073666384 , stdev= 102.12372098425168 Pecan Trees: avg= 322.8010876132931 , stdev= 178.50904048131653 Forest Trees: avg= 506.33531500288996 , stdev= 178.50904048131653 NIR band stats using masks Houses roof: avg= 7291.050868680545 , stdev= 985.2079469706347 Water: avg= 70.11727349703641 , stdev= 14.895863686638762 Pecan Trees: avg= 4829.812930513595 , stdev= 588.8033590079334 Forest Trees: avg= 4845.084964082239 , stdev= 588.8033590079334 SWIR1 band stats using masks Houses roof: avg= 6484.887462826733 , stdev= 907.4672793325856 Water: avg= 4.6685012701100765 , stdev= 7.364797051045533 Pecan Trees: avg= 1905.8157099697885 , stdev= 498.6547040104472 Forest Trees: avg= 1957.9827429609445 , stdev= 498.6547040104472 SWIR2 band stats using masks Houses roof: avg= 4933.191579276882 , stdev= 712.961764485916 Water: avg= 4.484335309060119 , stdev= 6.789876423663231 Pecan Trees: avg= 578.5833232628398 , stdev= 284.8599892051938 Forest Trees: avg= 655.9995871521758 , stdev= 284.8599892051938
#Using the predifined mask -
# Mask values
# 1 Water
# 2 White roofs
# 3 Walnut trees
# 4 Forest
# 5 Grass
# 6 Road
# 0 everything else
Objects_Label = ['Water','White Roofs','Pecan Trees','Forest','Grass','Road','party (everything else)']
Objects_Color = ['blue','Black','brown','darkgreen','lime','grey','red']
Avg_Spectra=np.zeros((7,426))
for band in range(0,nbands) :
Data2D=DataCube[:,:,band]
Avg_Spectra[0,band]=np.average(Data2D[ MaskData2 ==1 ])
for band in range(0,nbands) :
Data2D=DataCube[:,:,band]
Avg_Spectra[1,band]=np.average(Data2D[ MaskData2 ==2 ])
for band in range(0,nbands) :
Data2D=DataCube[:,:,band]
Avg_Spectra[2,band]=np.average(Data2D[ MaskData2 ==3 ])
for band in range(0,nbands) :
Data2D=DataCube[:,:,band]
Avg_Spectra[3,band]=np.average(Data2D[ MaskData2 ==4 ])
for band in range(0,nbands) :
Data2D=DataCube[:,:,band]
Avg_Spectra[4,band]=np.average(Data2D[ MaskData2 ==5 ])
for band in range(0,nbands) :
Data2D=DataCube[:,:,band]
Avg_Spectra[5,band]=np.average(Data2D[ MaskData2 ==6 ])
for band in range(0,nbands) :
Data2D=DataCube[:,:,band]
Avg_Spectra[6,band]=np.average(Data2D[ MaskData2 ==0 ])
plt.figure(figsize=(12,6))
plt.title('Features Spectral Signatures',fontsize=20)
plt.xlabel('Wavelength (nm)',fontsize=12)
plt.ylabel('Reflectance ',fontsize=12)
for i,type in enumerate(Objects_Label):
Spectra=Avg_Spectra[i]
plt.plot(Xvalues,Spectra,color=Objects_Color[i],label=Objects_Label[i])
plt.grid()
plt.legend(loc='upper left')
<matplotlib.legend.Legend at 0x7f9856796b80>
Going to generate a mask for only schrub and grassland based on this publication of NDVI value ranges https://www.researchgate.net/figure/Suitable-normalized-difference-vegetation-index-NDVI-ranges-identified-for-the-land_tbl1_330284135
# Mask Schrub and grassland (using NDVI band) 0.18-0.27
MaskSL= (NDVI > 0.18) & (NDVI < 0.27)
fig, figarr = plt.subplots(1, 2, figsize=(12,6))
fig.tight_layout()
figarr[0].set_title("Schrub and Grassland Mask using NDVI 0.18-0.27 range")
figarr[0].axis('off')
figarr[0].imshow(MaskSL)
# Display RGB Image
RGBImage=vip.Image_getRGB(DataRED,DataGREEN,DataBLUE,8000)
# Display RGB True color Image
figarr[1].axis('off')
figarr[1].set_title("NEON RGB - This a composite image made up of pieces")
figarr[1].imshow(RGBImage)
<matplotlib.image.AxesImage at 0x7f9856f5bc70>
Sort of did this above but now I'm going to see if I can do this just for urban sprawl (same source as above) with NDVI threshold of 0.015–0.14
# Mask urban sprawl (using NDVI band) 0.015–0.14
mmin=0.015
mmax=0.14
Name = 'Urban Sprawl'
Mask= (NDVI > 0.015) & (NDVI < 0.14)
fig, figarr = plt.subplots(1, 2, figsize=(12,6))
fig.tight_layout()
figarr[0].set_title(Name+' mask using NDVI '+str(mmin)+'-'+str(mmax)+' range')
figarr[0].axis('off')
figarr[0].imshow(Mask)
# Display RGB Image
RGBImage=vip.Image_getRGB(DataRED,DataGREEN,DataBLUE,8000)
# Display RGB True color Image
figarr[1].axis('off')
figarr[1].set_title("NEON RGB - This a composite image made up of pieces")
figarr[1].imshow(RGBImage)
<matplotlib.image.AxesImage at 0x7f985772e340>
#Using the predifined mask -
# Mask values
# 1 Water
# 2 White roofs
# 3 Walnut trees
# 4 Forest
# 5 Grass
# 6 Road
# 0 everything else
CustomMask= (NDVI > 0.015) & (NDVI < 0.14)
Objects_Label = ['White Roofs','Road','Custom Urban Mask']
Objects_Color = ['Black','grey','red']
Avg_Spectra=np.zeros((3,426))
for band in range(0,nbands) :
Data2D=DataCube[:,:,band]
Avg_Spectra[0,band]=np.average(Data2D[ MaskData2 ==2 ])
for band in range(0,nbands) :
Data2D=DataCube[:,:,band]
Avg_Spectra[1,band]=np.average(Data2D[ MaskData2 ==6 ])
for band in range(0,nbands) :
Data2D=DataCube[:,:,band]
Avg_Spectra[2,band]=np.average(Data2D[CustomMask])
plt.figure(figsize=(12,6))
plt.title('Features Spectral Signatures',fontsize=20)
plt.xlabel('Wavelength (nm)',fontsize=12)
plt.ylabel('Reflectance ',fontsize=12)
for i,type in enumerate(Objects_Label):
Spectra=Avg_Spectra[i]
plt.plot(Xvalues,Spectra,color=Objects_Color[i],label=Objects_Label[i])
plt.grid()
plt.legend(loc='upper left')
<matplotlib.legend.Legend at 0x7f98576c0ac0>
The spectral signatures of my custom mask versus the other urban values is in between the two other spectral signatures which is exactly what I would expect since it is averaging across all urban land. Also, I would expect its reflectance to be much less than the white roofs themselves but since it includes the white roofs, I would expect it to be more than roads. Overall, I'm pretty pleased with this mask.
#display a message the program ended
print("program ended.")
program ended.